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ABSTRACT 

Recent numerical relativity simulations have shown that the emission of gravitational waves during the 
merger of two supermassive black holes (SMBHs) delivers a kick to the final hole, with a magnitude as large 
as 4000km s"'. We study the motion of SMBHs ejected from galaxy cores by such kicks and the effects on 
the stellar distribution using high-accuracy direct A^-body simulations. Following the kick, the motion of the 
SMBH exhibits three distinct phases. (1) The SMBH oscillates with decreasing amplitude, losing energy via 
dynamical friction each time it passes through the core. Chandrasekhar's theory accurately reproduces the 
motion of the SMBH in this regime if 2 < In A < 3 and if the changing core density is taken into account. (2) 
When the amplitude of the motion has fallen to roughly the core radius, the SMBH and core begin to exhibit 
oscillations about their common center of mass. These oscillations decay with a time constant that is at least 10 
times longer than would be predicted by naive application of the dynamical friction formula. During this phase, 
the SMBH is typically displaced from the peak of stellar density by roughly the core radius. (3) Eventually, 
the SMBH reaches thermal equilibrium with the stars. We use straightforward scaling arguments to estimate 
the time for the SMBH's oscillations to damp to the Brownian level in real galaxies and infer times as long 
as ~ 1 Gyr in the brightest galaxies. The longevity of the oscillations makes this mechanism competitive with 
others that have been proposed to explain double or offset nuclei. Ejection of SMBHs also results in a lowered 
density of stars near the galaxy center; mass deficits as large as five times the SMBH mass are produced for 
kick velocities near the escape velocity. We compare the A^-body density profiles with luminosity profiles of 
early-type galaxies in Virgo and show that even the largest observed cores can be reproduced by the kicks, 
without the need to postulate "hypermassive" binary SMBHs. Implications for displaced AGNs and helical 
radio structures are discussed. 
Subject headings: galaxies: nuclei - stellar dynamics 



L INTRODUCTION 

The recent breakthroughs in numerical relativity dPretoriusI 
I2005t ICampanelh etaP 120061 iBaker et alJ r2006ah have al- 
lowed a number of groups to evolve binary black holes (BHs) 
to full coalescence. The final inspiral is driven by emission of 
gravitational waves, and in typical (asymmetric) inspirals, a 
net impulse is impart ed to the system due to anisotropic emis- 
sion of the waves jBekensteinl 119731; iFitchett & Detweile^ 
[THi [Favata et al."2004'). Early arguments that the magni- 
tude of the recoil velocity would be modest for non-spinning 
BHs jRedmount & Reesll9 89) were confirmed by the simula- 
tions, which found Vkick 5, 200kms~' in the absence of spins 
(iBaker et al.1 l2006bl; lOonzalez et alJ l2007bt iHerrmann etaP 
l2007h . The situation changed dramatically following the 
first simulations of "generic" binaries, i.e., binaries in which 
the individual BHs were spinning and in which the spins 
were allowed to have arbitrary orientations (ICampanelh et alJ 
l2007b) . Kicks as large as ^ 2000 kms~' have now been 
confirmed (Campanelli e t al.l l2007at iGonzalez et al.l l2007al ; 
iTichy & Marronetti 2007), and simple scaling arguments sug- 
gest that the maximum kick velocity would probably increase 
to ~ 4000 km s~' in the case of maximally-spinning holes 
dCampanelli et al.ll2007al) . The most propitious configuration 
for the kicks consists of an equal-mass binary in which the in- 
dividual spin vectors are oppositely aligned and oriented par- 
allel to the orbital plane. The kick amplitude also depends 
sensitively on the angle between the BH spi n vectors and their 
linear momenta shortly before the plunge (ICampanelh et alJ 
120073) . 

Electronic address: ' alessiag,merritt@astro. rit.edu | 



Galaxy escape velocities are < 3000 km s ^ dMerritt et alj 
|2004), which means that gravitational wave recoil can in prin- 
ciple displace coalescing supermassive black holes (SMBHs) 
arbitrarily far from galaxy centers, or even eject them com- 
pletely. The actual distribution of kick velocities is very 
uncertain, since it depends on the unknown distribution of 
binary mass ratios and spins, but most kicks are probably 
< 10 kms"'. A SMBH that is kicked with less than escape 
velocity will travel some maximum distance from the galaxy 
center after which its orbit decays due to dynamical friction; 
most of the energy loss takes place during passages through 
the galaxy center. Removal of the SMBH from the core has 
the effect of transferri ng kinetic energy to the stars and low- 
ering t he core densitv ((Redmoun & Reeslfl989l ; lMerritt et all 
2004; Boylan-Kolchin et alJl2004l) . This implies a more grad- 
ual return of the SMBH to a zero-velocity state than in a 
galaxy with fixed density. 

In fact, however, the SMBH is not expected to ever reach a 
state of zero kinetic energy. When its energy falls to a value 



1 , 1 , 



(1) 



with respect to the galaxy central potential, where m,^ and 
V* are a typical stellar mass and velocity respectively, ran- 
dom gravitational perturbations from stars act to accelerate 
the SMBH as often as they dec elerate it. This is the regime of 
gravitational Brownian motion ( lYounglll977l ; lBahcall & Woi3 
11976; Merritt et al. 2007). A natural definition of the "re- 
turn time" of a kicked SMBH is the time required for dy- 
namical friction to reduce the SMBH's mean kinetic energy 
to the Brownian value. Applying standard expressions for 
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the dynamical friction force leads one to the conclusion that 
this would occur in a relatively short time, of order a few or- 
bital periods, after dynamical friction has returned the kicked 
SMBH to the core. 

The A^-body simulations presented here were designed to 
test these expectations by evaluating the return times of kicked 
SMBHs and by quantifying the induced changes in galaxy 
structure. These processes can not be studied accurately using 
classical dynamical friction theory since the SMBH substan- 
tially modifies the core as it recoils and falls back. Approx- 
imate A^-body schemes, e.g. tree or grid codes, are also not 
well suited to the problem since they can not robustly follow 
both the early (collisionless) and late (collisional) evolution 
of the SMBH. Large particle numbers are required in order to 
cleanly separate the collisional and collisionless regimes. 

These various requirements can currently be met only 
with parallel, direct A^-body codes running on special- 
purpose supercomputers. Ou r sim ulations use the 

0GRAPE integrator (iHarfst et al.1 l2007h as implemented 
on gravitySimulator, a 32-node supercomput er em- 
ployi ng GRAPE-6A accelerator boards ( Fukushi ge et al.l 
I2005h . 

Our findings are surprising in one important respect. After 
returning to the core, the kicked SMBH exhibits long-lived 
oscillations with amplitude comparable to the core radius'. 
These oscillations eventually decay but with a time constant 
that is at least an order of magnitude longer than would be 
predicted by a straightforward application of the dynamical 
friction equation. We demonstrate that the existence, ampli- 
tude and damping time of these oscillations are independent 
of the number of "star" particles used in the simulations, for 
A^ up to 2 X 10*'. The oscillations are si milar to those first re- 
ported by R Miller and collaborators ( Miller & Smithll 19921: 
liviiUeilfT996h in their pioneering A^-body studies of the central 
regions of galaxies. A number of other authors have reported 
low effective values of the dynamical friction force as it acts 
on massive objects that inspiral into constant-dens ity cores 
dBontekoel ligSS): iBertin et alj 120031: iRea d et al.H 2006^ or on 
rotati ng bars TWeinberg &'Katdl2002n Valenzue la & KlvpinI 
l2003h . Our use of a high-accuracy, direct-summation A^-body 
code combined with large particle numbers greatly reduces 
the possibility that our results are an artifact of the potential 
calculation scheme, an issue that has plagued the interpr eta- 
tion of similar results in the past (IZaritsky & Whitelll988h . 

§|2] describes the initial models and the A^-body algorithm. 
Evolution of the SMBH's orbit is described in detail in §[3] 
and the induced changes in galaxy structure are described in 
§21 where the A^-body models are compared to luminosity 
profiles of core galaxies. §|5]presents estimates of the SMBH 
return times in real galaxies, and §|6] discusses some of the 
observable consequences of the kicks. 

2. INITIAL MODELS AND NUMERICAL METHODS 

The light profiles of elliptical galaxies and the bulges of spi- 
ral galaxie s are general l y well desc ribed in terms of the Ser- 
sic model ^Sersidll963t I Sersi3ll968 ). which is a generaliza- 
tion of thejde Vaucouleurs] (ll948l [1959) law. The most lumi- 
nous elliptical galaxies depart systematically from the Sersic 
law near the center, where t hey s how evide nce for partially 
deple t ed stellar cores (Faberetal. 1997: Mi losavlievic et all 
120021: iGrahamI [2004; .Ferrarese et al.. ,2006.) . Formation of 

* A movie showing the oscillatio ns is av ailable at 
|http://ccrg . rit.edu/Research/Publications. php?pap er=0708.077r[ 



TABLE 1 

Parameters of the initial models. 



name 


n 


a 


'•b 


7 




Al 


4.0 


2.0 


0.014 


0.55 


1.0 X 10-3 


A2 


4.0 


2.0 


0.0095 


0.55 


1.0 X 10-3 


B 


4.0 


2.0 


0.027 


0.55 


3.0 X 10-3 



a binary SMBH following a galaxy "major merger" has 
been shown to produce cores of roughly the right magnitude 
(Milosavlievic & Merritt 2001: Merritt 2006), although some 
observed cores are too large to be easily explained by this 
model (a point we return to in detail below). 

As approximate representations of galaxies with binary- 
depleted cores, we adopt core-Sersic models ([Graham et alj 
2003) for our initial conditions. The space density profile of 
a galaxy that follows the co re-Sersic law in projectio n can be 
accurately approximated as ([Terzic & Grahamll2005[) 



with 

p =Pb2<''-^'/"(^^y/('""^'/"")"". (3) 

Equation ([Hi is a modifica tion of the Prugniel-Simien model 
( [Prugniel & Simied [T997[) . Here, is the effective (half- 
mass) radius of the projected galaxy: '"b is the break (core) 
radius: Ph is the space density at r= r^: and a regulates the 
sharpness of the transition from core to outer profile. The 
parameter n describes the curvatur e of the Sersic profile and 
b and p are fixed functions of n (Prugniel & Simienlll997[ : 
[Terzic & Graham 2005). Monte-Carlo initial conditions were 
generated using the scheme of Szell et al.l (1200 5). after includ- 
ing the gravitational potential of a central point particle repre- 
senting the SMBH. 

The parameters used for our initial models are listed in Ta- 
ble[T] The table also reports names for the different runs based 
on the adopted ratio of SMBH mass to galaxy mass and ini- 
tial core radius. Core radii were chosen so as to give initial 
mass deficits of roughly Mbh, as observed fo r the majority 
of luminous early-type galaxies (' Merritl[[2006 f). We note that 
7 = 0.5 is the shallowest power-law profile that is consistent 
with a non-negative, isotropic distribution of stellar velocities 
around the BH. 

The initial models were evolv e d usin g the (/iGRAPE 
numerical integrator ([Harfst et alj [20071) . This direct- 
summation code employs a fourth-order Hermite integrator 
with predictor-corrector scheme and hierarchical time steps. 
The MPI parallelization strategy is designed to minimize the 
amount of communication among different computing nodes 
and to make efficient use of the special-purpose GRAPE hard- 
ware. All the simulations presented in this work were per- 
formed on the 32-node cluster gravitySimulator at 
the Rochester Institute of Technology. Most of our simula- 
tions used A^ = 0.5 x 10^ equal-mass particles to represent the 
galaxy although some runs used larger A^. We set the ratio 
of BH mass to galaxy mass, Mbh/Ms^i, to be (1,3) x 10~^, 
typical for observed galaxies (Merritt & Ferrarese '2001*). For 
each model described in Table [U we chose eleven different 
values of the kick velocity Vkick in units of the central escape 

^ http://wiki.cs.rit.edu/bin/view/GRAPEcluster 
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speed Vesc: Vidck = (0.1,0.2, .., 1.1) x Vesc; the latter was com- 
puted numerically from the initial A^-body models. In order 
to guarantee energy conservation, we used a time-step accu- 
racy parameter r/ = 0.01. This ensures a relative energy error 
smaller than one part in 10^. An accuracy parameter twice as 
big would approximately halve the integration time but would 
result in a relative energy error of 10"^, which we do not con- 
sider acceptable for this study. A softening length e = 10"^ 
was assigned to both the stars and the BH. Such a small soft- 
ening length has been shown not to affect even the Br ownian 
motion of a massive particle in models like ours ( Merritt et alj 
lIOOTh . 

Throughout the paper we adopt units according to which the 
gravitational constant G, the effective radius 7?e in equation 
(|2|i, and the total galaxy mass Mgai are unity. The models can 
be scaled to physical units as follows; 



[T]: 



-1/2 



7.75 X lOV 

1/2 



101' M, 



-1/2 



3kpc 



3/2 



[V] = 



GM. 



gal 



= 378kms- 



gal 



lOiiM 







Ukpcj 



1/2 



(4) 



(5) 



(6) 



(7) 



3. THE BLACK HOLE MOTION 
3.1. General Remarks 

Figure [H (which ca n be compared with Fig. 1 of 
iMadau & Ouataerl2004l shows BH trajectories in models Al, 
A2 and B for Vkick/Vesc = (0.4,0.7,0.9, 1.1). For Vyck > Vesc 
the black hole escapes the galaxy on an unbound orbit. The 
maximum displacement of the BH (r^ax) is shown in Figure|2] 
The data from the simulations (points) are compared to theo- 
retical (dotted lines) and numerical (dashed lines) estimates of 
'"max in the absence of dynamical friction. The theoretical and 
numerical estimates are obtained from the initial A^-body data 
by assuming conservation of total energy for the BH: r^ax is 
the distance at which the gravitational potential of the system 
equals the initial total energy of the BH. For the theoretical so- 
lution we use the expression of the potential in a core-Sersic 
model (see equations 7 through 13 of Terzic & Graham 2005|) 
while for the numerical solution we compute the potential at 
different radii from the A^-body data. The two estimates are 
for practical purposes indistinguishable. 

Dynamical friction affects the maximum displacement of 
the BH only for moderately large kicks, where the data points 
appear systematically lower than the theoretical curves. Val- 
ues of Tmax larger than the expected turning points in the first 
orbit are due to the rapidly-expanding core. 

During the initial outward journey, dynamical friction does 
not strongly influence the motion of the BH, and the maxi- 
mum displacement is similar to that of an energy-conserving 
orbit. We note that a kick velocity larger than about 0.3 Vesc 
is necessary to bring the BH beyond the core. Due to the 
combined effect of the kick and dynamical friction, the BH 
displays a damped oscillatory motion. The number of radial 
oscillations increases with Vidck; for Vkick = 0.9 Vesc the BH ex- 
periences ~ 5 full radial oscillations before returning to the 
core. Almost all of the energy loss to dynamical friction takes 
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Fig. I. — BH trajectories in models Al, A2 and B, for Vkick/K:s 
(blue/lower), 0.7 (green), 0.9 (red) and 1.1 (black). 
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place during the short intervals that the BH passes through the 
core. This is shown in Figure |3] which plots the evolution of 
the BH specific energy E in Model Al with V^ck = 0.9Vesc, 
where 

T/2 ^ 

^ -E . (8) 



E = 



V(x,-X)2 + e2 
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Fig. 2. — Maximum displacement of the BH from tlie galaxy center. The 
data points show the results from the simulations while the lines are estimates 
in the absence of dynamical friction. The dashed lines represent numerical 
estimates from the computation of the potential of the W-body system at time 
^ = while the dotted lines represent theoretical estimates from the analytic 
expression of the potential in a core-Sersic model. 
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Fig. 3. — Upper panel: Specific energy of the BH particle versus time in 
Model Al with Vi;i^.ii = 0.9V'esc. Almost all of the energy loss occurs during 
passages through the core. Lower panel: Mean density in a sphere of radius 
0.05 centered on the point of maximum density in the core of the galaxy 
(excluding the BH). 

and the summation is over the "star" particles^. The energy 
lost during the initial emergence from the core appears to be 
less than during subsequent passages, suggesting that dynam- 

^ Unless otherwise noted, upper-case variables X and V refer to the BH 
particle while lower-case symbols are reserved for the star pailicles. 



Fig. 4. — Stellar mass bound to the BH in the initial models. 

ical friction requires a finite time to "turn on" after the kick. 
During the first few oscillations, the BH's motion remains es- 
sentially rectilinear, but eventually the Y- and Z-components 
of the motion become important due to non-sphericities in the 
galaxy potential and also to perturbations from stars. At late 
times, the BH's motion is essentially random, similar to that 
of a Brownian particle in a fluid. Figure [3] also shows the 
mean density in a sphere of fixed radius whose center is lo- 
cated at the estimated density peak (computed via the algo- 
rithm described in § 13.21 ). The core density decreases rapidly 
following the initial ejection, then more gradually as the BH 
returns again and again to the core, losing energy to the stars 
each time. 

Figure m shows the mass in stars bound to the BH at f = 0. 
The bound mass was computed by counting all the stars, 
within the influence radius r/,, which formed a bound two- 
body system with the BH particle. The influence radius was 
defined as the radius containing a mass in stars equal to twice 
Mbh- The bound mass decreases steeply with Vkick, as noted 
in ea rlier studies (Merritt et al. 2004; Bovlan-Kolchin et al] 
|2004|) . and is ignorable for Viick ^ 0.6yesc- 

In all cases where the kick velocity was large enough to re- 
move the BH completely from the core (i.e. Vkick ^ O.SVesc), 
we observed three distinct regimes of the motion. In Phase 
I, the BH's motion is well predicted by Chandrasekhar's dy- 
namical friction theory, after taking into account the chang- 
ing size of the galaxy core where most of the friction occurs. 
This is the phase illustrated in Figure [TJ in Figure [3j Phase I 
extends until t w 20. Phase II begins roughly when the ampli- 
tude of the BH's motion had decayed to the size of the core. 
In this phase, the energy of the BH's orbit continues to decay 
but with a much longer time constant than predicted by Chan- 
drasekhar's formula. The BH and the core oscillate about their 
common center of mass in this regime. In Phase III, the BH's 
energy has dropped to the thermal level. Phase II is generally 
longer than Phase I, and this would presumably be even more 
true in real galaxies since the amplitude of thermal oscilla- 
tions is much lower than in our simulations implying a longer 
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time to reach the Brownian regime. We discuss these three 
regimes in detail below. 

3.2. Phase I 

The extent of Phase I is clearly indicated in the plots of BH 
energy vs. time (e.g. Fig.|3]i: a distinct "knee" appears in the 
E(t) curves marking the end of this phase. Values of T/, the 
elapsed time from the kick until the end of Phase I, are given 
in Table |2] 

We compared the evolution of the BH's motion in Phase 

1 with the predictions ofChandrasekhar's dynamical friction 
theory (ChandrasekJi^ [19431). Such comparisons are prob- 
lematic since much of the energy exchange between BH and 
stars occurs during passages through the galaxy's core, and 
the core density changes significantly with time due to the 
BH's motion. We dealt with this problem by breaking the 
BH's motion into segments, each containing one passage 
through the center, and assuming that the galaxy's density re- 
mained constant during each segment. 

Chandrasekhar (1943) derived his expression for the dy- 
namical friction acceleration F^f assuming an infinite, homo- 
geneous and unchanging background of perturbers (stars). In 
the limit that the mass of the heavy object greatly exceeds the 
masses of the stars, the acceleration is predicted to be 

Fdf w -iTTG^pMYia \n{\ + K^)V-^N{< y,r), (9) 

where p(r) is the mass density of stars at the BH's position, 
(1 /2)ln(l + A^) is the Coulomb logarithm, V is the BH's in- 
stantaneous velocity, and A^(< y,r) is the fraction of stars at 
r that are moving (in the frame of the galaxy) with velocities 
less than V . 

Some care must be taken in the definition of the Coulomb 
logarithm. One commonly writes 

ln(l + A2) «21nA«21n(p„ax//?min) (10) 

where pmh, and p^ax are the minimum and maximum effective 
impact parameters of the stars that contribute to the frictional 
force, and P max ^ Pmin- However, Pmin depends on the field- 
star velocity (IWhitelll949t lMerrittll200lb and pmax is likewise 
ill-defined since a realistic stellar system is inhomogeneous 
and has no outer boundary. 

Numerous A^-body simulations have been carried out to 
evaluate Chandrasekhar' s formula in the case of a mas- 
sive particle inspiraling toward the center of a galaxy 
jWhite 1983; Bontekoe & van Albada 1987; Bontekoe 1988; 
I^inberg|ll989t ICora et al.lll997l: iBertin et all 120031) . Early 
work was typically based on approximate A-body schemes 
and the results were of ten discrepant from study to study 
dZaritskv & White! 1 198 8*). These differences appear to have 
been resolved in the last few years through the use of direct- 
summation codes ( S pinnato et al.ll2063l: iMerr itt 2006), which 
consistently find 4 < InA < 6 for inspiral of massive point 
particles, on ckcular or near-circular orbits, into the centers 
of galaxies with steeply-rising density profiles. Fewer exper- 
iments have been done w ith highly eccentric orbits, although 
lJust & Peiiarrubial (l2005h . using an approximate method, find 

2 < In A < 3 for orbits with moderate eccentricities. 

In general, we expect the effective value of InA to be 
smaller for radial orbits than for circular motion. The dy- 
namical friction force arises from a polarization of the stellar 
density which produces an over-dense region, or wake, be- 
hind the massive object (Mulder 1983i) . A finite time, of order 
a galaxy crossing time, is presumably required for this wake 



to be set up. In the case of a gradually-decaying circular or- 
bit, the galaxy is able to reach a quasi-steady state after a few 
orbits of the massive object. In our case, the position and 
velocity of the BH are changing dramatically over one cross- 
ing time, so that the wake never has a chance to establish its 
steady-state amplitude; indeed just after apocenter passages, 
the over-dense region can be seen to lie in front of the BH. 

In order to determine the effective value of In A in the A- 
body integrations, we computed BH trajectories using Chan- 
drasekhar's formula (equation|9]) with various values of the 
ln(l+A^) term (henceforth written simply as 2 In A) and 
compared them with the A-body trajectories. The following 
procedure was followed. 

1. The density center of the galaxy moves slightly with 
respect to the origin of the coordinates due to transfer of mo- 
mentum from the kicked BH to the galaxy. In order to accu- 
rately determine the distance of the BH from the galaxy center 
as a function of time, we recorded full snapshots of the par- 
ticle positions at frequent intervals, then used the Casertano- 
Hut (1985) algorithm to find the density center of the stars 
in each snapshot. A smoothing spline was fit through the 
measured positions to give a continuous estimate of the center 
displacement as a function of time, and this displacement was 
subtracted from the BH positions. (The instantaneous velocity 
of the density center was ignored, which is a good approxima- 
tion at least until the end of Phase I.) The resulting correction 
was at most ~ 0.02; at late times the displacement reached a 
constant value since the center-of-mass velocity of the system 
was zero by construction. 

2. In order to apply Chandrasekhar's formula we needed 
to specify the galaxy model. The galaxy's mass distribu- 
tion changes with time due to the BH's motion; most of 
this change takes place in the core just after the BH passes 
through. We therefore fixed all the parameters in equation (|2]i 
except for the core radius r/,. We determined the effective 
value of rh at the discrete times when the BH passed through 
the galaxy center by assuming a flat core (7 = 0) and finding 
the value of r/, such that the mass contained within rh accord- 
ing to equation (|2]l, with a = 2, was the same as the mass in 
the A-body model in a sphere of radius rh centered on the BH. 
This procedure was always found to yield a unique rh and ac- 
curately recovered the known value of rh in the initial models. 

3. BH trajectories were then computed in a piece wise 
fashion using Chandrasekhar's formula, starting from one ex- 
tremum in the BH displacement and continuing until the next 
extremum, using the value of rh corresponding to the central 
passage lying between the two extrema. This w as repeated 
for sev eral values of InA. We used equation (5) of lSzell et alJ 
( l2005h to compute A(< V, r) in equation ^ from the assumed 
P(r). 

Figure |5] shows the results for Model Al with Vkick/Vesc = 
0.7 and Model B with Vkick/V'esc = 0.8. During each inward 
leg of the trajectory, the dynamical friction force hardly af- 
fects the motion; only when passing through the dense center 
is the motion significantly non-ballistic. (This could be seen 
already in Figures|2] and [3]) The best-fit value of InA was 
found to lie in the range 2 < InA < 3, and for such values, 
Chandrasekhar's formula did a good job of reproducing the 
motion. We found no evidence of a systematic change in the 
effective value of In A from one time interval to the next. 

3.3. Phase HI 

The BH trajectories in Figure |5] are displayed until the am- 
plitude of the oscillations has decayed down to roughly the 
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Fig. 5. — Comparison between BH trajectories computed via the W-body 
integrations (open circles) and via Chandrasekhar's formula (9) (lines). The 
/V-body models were Al (Mbh = 0.001) with Vu^k = 0.7 Vcsc (a) and B (Mbh = 
0.003) with Vyi-k = 0.8 Vesc (b). Theoretical trajectories were computed in a 
piecewise manner, starting from extrema in the BH's trajectory (vertical solid 
lines) and continuing until the next extremum; the core radius of the galaxy 
model was adjusted as described in the text to give the same core density as 
in the W-body model at the time when the BH passed through the center. 
Horizontal dashed lines show the adopted values of ri,. Line colors/styles 
correspond to different values of In A: 1 (blue/solid), 2 (magenta/dashed), 3 
(red/dash-dotted), 4 (black/dotted). 

core radius. As discussed above, the BH's motion is well pre- 
dicted by Chandrasekhar's dynamical friction formula in this 
regime. Shortly after returning to the core, however, the BH's 
motion was found to depart strikingly from the predictions of 
Chandrasekhar's formula. A detailed discussion of the mo- 
tion in "Phase 11" is presented below. Before doing so, we 
consider the motion of the BH at still later times, "Phase III," 
when it has reached thermal equilibrium with the stars. 

Figure |6] shows the squared velocity of the BH, = + 
Vy + y/, over the full integration interval, for kick veloci- 
ties Vkick > O.SVesc in Model B. For V^ick > O-^V^.c the BH 
moves substantially beyond the core during its first oscillation 
(Fig.|2]i. At late times, the motion of the BH in each of these 
integrations appears to be stochastic (i.e. non-quasi-periodic) 
but with roughly constant amplitude. 

The dashed (blue) lines in this figure show (V^), the mean 
square velocity of the BH averaged over Phase III. (The pre- 
cise definition of the start of Phase III is given below.) Also 
shown (dotted red lines) are estimates of the expected value 
of ( V^) for the BH once it reaches statistical equilibrium with 
the stars. The latter velocity, VQ^.^^^, was computed using 
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Fig. 6. — Squared BH velocity in seven W-body integrations of Model B. 
For Vkick 0.4Vesc the BH moves completely out of the core before falling 
back. Ticked, horizontal lines demarcate Phase II. Blue (dashed) lines show 
(V-) during Phase III, and red (dotted) lines show the mean square velocity 
predicted by equation jilt , which assumes that the BH particle has reached 
thermal equiUbrium with the stars in its vicinity. 



Equation ( fTTT i equates the kinetic energy of the BH with the 
mean kinetic energy of a single star in the core. The quantity 
(T is defined as the ID velocity dispersion of stars within a 
sphere of radius K x r/, centered on the BH, with r/j the BH's 
influence radius (the radius containing a mass i n stars equal 
to twice Mbh) and K a constant of order unity. iMerritt et alj 
((2007) used A^-body simulations to evaluate K for massive 
particles at the centers of galaxies with power-law nuclear 
density profiles, p ^ r~'' . They found that K increases slowly 
with decreasing j, to K ^ 0.8 when 7 = 0.5. We set A" = 1 
when computing Vbiowii in Figure |6j the agreement with the 
measured values is quite good, confirming that the BH be- 
haves as a Brownian particle in Phase III. 

Figure |7] shows the rms amplitude of the BH's motion aver- 
aged over Phase III. Since the density center of the galaxy 
drifts, as described above, smoothing splines were first fit 
to the X(t) values for the BH and the rms deviations were 
computed with respect to the smoothed trajectories. Figure|7] 
shows a general trend of increasing R^ms with decreasing core 
density, as expected if the motion in this regime obeys the 
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Fig. 7. — RMS amplitude of the BH oscillations in the Brownian 
regime. Phase III, for models Al (black/circles), A2 (blue/squares) and B 
(red/triangles). 



virial theorem. 



(y2) « -TrGpciR^ 



(12) 



This relation (cf. iBahcall & Woig iT976^ assumes a constant- 
density core, ignores the back-reaction of the BH's motion 
on the stars, and ignores any coupling between random grav- 
itational perturbations from the stars and the quasi-periodic 
motion of the BH in the smooth potential of the core. Nev- 
ertheless, equation ( fT2] l was found to reproduce the measured 
7?rms values in Figure [7] quite well if pc was defined as the 
mean density of stars within r/,. Fluctuations in /Jims about the 
mean relation in Figure |7] appear to be due primarily to fluc- 
tuations in Vtms and would presumably be smaller if the /?rms 
values were averaged over longer time intervals. The near 
agreement between the /?nns values for the runs with small 
and large Mbh is a consequence of the larger core size / lower 
core density in runs with larger Mbh, which compensates for 
the lower (V^) cx M'^l^. 

We note here that the amplitude of the BH's Brownian mo- 
tion is always a factor 10 or more smaller than the final core 
radii of the models (Table |3]l. This implies that the motion of 
the BH when it first returns to the core - at the start of Phase 
II - should not be appreciably affected by discreteness effects, 
i.e. by perturbations from individual stars. This conclusion is 
confirmed below. 

We note also that the amplitude of Brownian oscillations 
of BHs in real galaxies (expressed as a fraction of the galaxy 
effective radius, say) would be smaller than in our models by 
the factor ^ y^(M^^Jm^JjN , i.e. ^ 50 for Mgai = lO^M© and 

500 for Mgai =10' 'M0 . The time required for a BH to reach 
these lower kinetic energies would also presumably be longer 
than in our simulations, as discussed in more detail below. 

3.4. Phase II 



As noted above, the motion of the BH after returning to the 
core, and before reaching the Brownian regime, is not well 
described by Chandrasekhar's formula. Here we consider the 
motion in this regime ("Phase 11"). 

Figure |6]reveals the following qualitative features. 

1 . The motion in Phase II is essentially oscillatory, with a 
period similar to that at the end of Phase I, i.e. roughly equal 
to the period of oscillation of a test particle moving in the 
stellar core. 

2. There is evidence of additional frequencies affecting the 
BH's motion. For instance, the amplitude of the oscillations 
sometimes appears to increase temporarily over several peri- 
ods in a manner suggestive of beats. 

3. Averaged over many periods, the mean amplitude of 
the oscillations decays, but with a time constant that is much 
longer than observed toward the end of Phase I. 

4. Near the end of Phase II, the motion becomes increas- 
ingly stochastic, presumably due to perturbations from indi- 
vidual stars. Eventually the BH rms velocity falls to the Brow- 
nian (thermal) level marking the start of Phase III. 

5. Phase II always begins roughly when the stellar mass in- 
terior to the BH's orbit is equal to Mbh. When V^ck S 0.3Vesc, 
the BH never escapes the core, and its motion appears to tran- 
sition directly from Phase I to Phase III. 

Based on Figure|6] the elapsed time in Phase II can be sub- 
stantially longer than the time spent in Phase I. Understanding 
the character of the motion in this regime is therefore crucial 
for predicting the expected displacement of a supermassive 
BH in a real galaxy following a kick. 

We begin by considering a simple model for damped oscil- 
lations of a massive particle in a constant-density core. While 
this model will fail to quantitatively reproduce the motion in 
Phase II, it provides a useful framework for discussing what 
is observed in the simulations. 

In the absence of dynamical friction, and neglecting the in- 
fluence of the massive particle's presence on core structure, 
the motion of the massive particle is simple harmonic oscilla- 
tion with frequency uj^ = ^ {Air /y)Gpc', Pc is the core density, 
assumed constant within a radius re- To this motion we add 
the acceleration due to dynamical friction. If the velocity dis- 
tribution of the stars that produce the friction is Maxwellian 
with ID velocity dispersion a^, and if the BH's velocity sat- 
isfies V <C cTc, the resulting equation of motion in any coordi- 
nate Xi is 

(13) 



where 



(14) 



TT G^pcMfinlnA 

is the dynamical friction damping time (lMerrittlll985l) . The 
condition for underdamped oscillations is 2uJcTdf > 1, where 



V6 



2 GVZpy^MBHlnA 



(15) 



(16) 



Mbh In A' 

with Mc = {4/3)TTpcrl the cor e mass; the second relation 
uses the "core-fitting" formula of lRood et al.l (Il972h . 



a^,=F^^Gpy, 



2 for our models. Thus 

2wc7df ~ 4 



Mr 



Mbh In A 



(17) 
(18) 
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Fig. 8. — Evolution of the BH kinetic energy in a series of integrations of 
model B with various N, and Vyck = 0.6Vesc. Top panel: Squared velocity of 
the BH versus time. Dashed lines at the right show the predicted values of 
in the Brownian regime (eq. llU . Bottom panel: Binned values of V- in 
Phases II and III. Dotted lines are least-squares fits to the binned data. These 
fits are plotted until the time at which they intersect the Brownian V^; these 
times are marked by the vertical solid lines. The latter are found to be spaced 
with roughly constant separation indicating that the time required for the BH 
to reach thermal equilibrium with the stars increases roughly as \nN. 

In our simulations (and in real galaxies), the right hand side 
of this expres sion is > 1, since core masses are ~ a few Mbh 
(lMerritll2006l) and 2 < In A < 3 (§|l2ll. It follows that the mo- 
tion of the BH should be under-damped, though not far from 
critically damped, after it re-enters the core. The solutions to 
equation (fTsT i in the under-damped regime are 



-'/27ir 



sin {ujct + (pi) ■ 



(19) 



Writing Q = luj^Tdf > 1 and 7^. = 2tt/uJc, the energy decay 
time is predicted to be T^f = {Q/ATi)Tc, i.e. shorter than the 
orbital period. Such short decay times are in fact observed 
near the end of Phase I (Figure |5]l. 

However, Figure|6]shows that this is not the case in Phase II: 
the mean damping time is substantially longer than an orbital 
period. The abrupt decrease in the energy dissipation rate at 
the start of Phase II can also be seen in Figure|3la). 

A possible explanation for the slower damping in Phase 
II is discreteness effects: perturbations from individual stars, 
some of which act to accelerate the BH, become increasingly 
competitive with mean-field effects (including dynamical fric- 
tion) as the BH moves more slowly. Indeed, in the Brownian 
regime (Phase III), the accelerating perturbations are equally 
as strong, in a time-averaged sense, as dynamical friction. 



While the ampUtude of the BH oscillations at the onset of 
Phase II is always much greater than the Brownian ampli- 
tude in these simulations (cf. Fig. |7] and the accompanying 
discussion), it is still conceivable that discreteness effects are 
responsible for the anomalously slow decay of the BH's orbit 
at this time. 

To securely rule out this possibility, we repeated the inte- 
gration of model B with Vyck = 0.6Vesc, increasing up to 
N = 2x 10^. Figure [8] shows the results. The slowly-damped 
oscillations in Phase II are clearly not an artifact of a too-small 
N. In all cases, for instance, the fifth extremum in (which 
occurs at f « 1 .42) is comparable or greater in amplitude to the 
fourth extremum (at f « 1.16), rather than being much lower 
in amplitude as would be expected from the above analysis or 
from Figure |5] We also carried out a number of tests varying 
the integration time-step parameter rj; again, no systematic 
dependence of the evolution in Phase II on this parameter was 
observed. 

Particularly striking in Figure [8] is the accurately exponen- 
tial decay of the BH's kinetic energy throughout Phase II; this 
is clearest in the simulation with largest N, where the expo- 
nential damping continues over two decades in energy. We 
note again that an exponentially decaying energy is predicted 
by the simple model just presented, but the model predicts a 
much shorter time constant than what is observed in the A- 
body simulations. 

Figure |9] suggests why Chandrasekhar's (1943) formula 
might break down in Phase II. The approximation of a station- 
ary galaxy is strongly violated in this regime. The galaxy's 
density center oscillates with opposite phase to the BH, and 
with roughly the same frequency and amplitude. This is 
consistent with the observation that Phase II always begins 
roughly when the mass in stars inside the BH's orbit is similar 
to Mbh. Evidently, in this regime, the BH and the core oscil- 
late about their common center of mass as a two-body sys- 
tem. Chandrasekhar's derivation, which assumed a body on 
a linear trajectory through an infinite homogeneous medium, 
is unlikely to apply to oscillations like those in Figure |9] 
since the BH is periodically accelerated, then decelerated, 
by the density peak. The rate at which such oscillations de- 
cay is known to be sen sitively dependent on resonant interac- 
tions dJ remaine & Weinberg 1984) an d can be arbitrari ly low 
(lLouis& G erhard 1988: Sridhar 1989t ISridhar & Nitvanandal 
1989; Mine au et al.lll99()h . although we are not aware of any 
theoretical treatment that is directly applicable to oscillations 
like those in Figure |9] 

Ours is not the first A-body study to observe persistent os- 
cill ations of massive obje cts a t the centers of A-body mod- 
els. iMiUer & SmIfliI (Il992h andlMiUeiJ (Il996h reported a series 
of A-body integrations, using a grid-based code, of a disk at 
the center of an axisymmetric galaxy model. They observed 
what appeared to be over-stable oscillations of particles ini- 
tially at rest near the center of the disk; the oscillation fre- 
quency was roughly y/ (4Tr/3)Gpc and the maximum ampli- 
tude was roughly the size of the core. All of these features are 
characteristic of the oscillations that we observe in Phase II. 
iMiller & Smithl(ll992i) also reported "a couple of experiments 
in which a massive object was put into orbit within a galaxy 
model," presumably near the center, and observed "residual 
oscillations" with amplitude roughly equal to the radius at 
which the enclosed mass wa s equal to the o bject's mass, again 
similar to what we observe. iMiller & SmithI (11992) briefly de- 
scribe a model for the oscillations, in which periodic motion 
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Fig. 9. — Core-BH oscillations in Phase II. This is the W = 2 X lO'' integra- 
tion of Model B shown as the filled (red) circles in Fig.[8] Contours are sep- 
arated by 0.034 in log[Q of the projected density. Filled circles mark the BH 
and crosses mark the approximate location of the (projected) stellar density 
maximum. Times are f = 2.1875,2.21875,2.25, ...2.46875, increasing from 
upper left to lower right. The elapsed time in this figure spans approximately 
1/2 oscillation period of the BH. 



of the core as a whole, at roughly the same frequency as core 
internal frequencies, drives the oscillations. 

A number of other A^-body studies have noted a decrease in 
the effective value of In A once a massive object has spiraled 
into a constant-density core. Typically the ob served decrease 
is modest, a factor^ 2-3 or so (.Bp ntekoe & van Albadal 
IT987t lBontekodlT988t IWeinberdUgsgTlC ora et al.l ll997h . al- 
though one recent studv dRead et al.l 2006) found a nearly 
complete disappearance of dynamical friction after the in- 
falling particle reached the core. Read et al. proposed that 
the apparent vanishing of the dynamical friction force in their 
simulations could be explained by the degeneracy of orbital 
frequencies in the harmonic-oscillator potential correspond- 
ing to a precisely flat, central density profile. In such models. 
Read et al. found that the disappearance of the dynamical fric- 
tion force was critically dependent on whether the plane de- 
fined by the inspiralling particle's orbit remained fixed; pre- 
cession, induced e.g. by finite-A^ perturbations, caused dy- 
namical friction to turn on again, though at a rate much slower 
than expected from Chandrasekhar friction. While Read et al. 
did consider the effect of varying the initial log-slope of the 
background distribution, their models were always spherical. 
The cores in our models are not precisely flat nor are our mod- 
els precisely spherical (once the BH particle has been ejected) 
and these differences (coupled with the fact that the gravita- 
tional potential of the core is highly oscillatory in Phase II) 
may explain why we do not observe the dramatic stalling re- 
ported by Read et al. In any case, the apparent lack of an 
A^-dependence in our simulations (Figure |8]i suggests that the 
critical difference between our results and those of Read et al. 
is not particle number 

The Phase II oscillations were clearly visible in every in- 
tegration with Vkick > 0.4Vesc- For Viick = O.SVesc there were 
hints of a delayed return to the Brownian regime in some of 
the integrations (e.g. Fig. |6]) but not to the extent that we were 
able to estimate damping times. We could not detect the Phase 
II oscillations at all for Vkick < 0.2V'esc; in these integrations, 
the BH kinetic energy appears to drop very rapidly after the 
return to the core, more or less as expected based on the an- 
alytic model presented above or by an extrapolation of the 
behavior in Phase I. In any case, we assume in the remain- 
der of this paper that the Phase II oscillations are absent when 
Vkick < 0.3Vesc- Integrations with much larger might modify 
this conclusion. 

The occasional increase in the amplitude of the Phase II 
oscillations, which is seen in virtuall y all the integrat ions, is 
suggestive of a dynamical instability (iTremainell2005h . How- 
ever an instability would presumably act even in the case of 
small kicks, while as noted above. Phase II oscillations appear 
to be absent for Vkick O.SVesc- We speculate that the BH must 
be kicked completely out of the core in order for the BH-core 
oscillations to be excited, as suggested by Mill er & SmithI 
[[992 ). The roughly sinusoidal variations in the envelope of 
V^{t), with a much lower frequency than ojc, could naturally 
be explained in terms of beating, e.g. between the frequency 
of motion in the core and the frequency at which the core itself 
oscillates in the galactic potential. 

Figures |6] and [8] suggest that the core-BH oscillations in 
Phase II decay roughly as an exponential in time, at least when 
viewed through a window of several orbital periods or longer. 
We investigated a number of ways to quantify the time con- 
stant r associated with the energy damping: 

1 . Plots of BH energy versus time (equation [8] Fig. O 
were found not to be very useful in this regard since the total 
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Fig. 10. — Energy-decay time constants r for the BH in Phase II, for 
models Al (black/circles), A2 (blue/squares) andB (red/triangles). 

energy is dominated by the potential energy which exhibits 
fairly large fluctuations from time step to time step. 

2. In a constant-density core, the unperturbed motion is 
simple harmonic oscillation with frequency oj and energy 

1 ^ 

We determined the dominant frequency of the BH's motion 
in Phase 11 by carrying out discrete Fourier transforms of the 
comple x functio ns Xi{t) + iViit) and constructing power spec- 
tra (e.g. lLaskaii ri990). Least-squares fits to In/ssHO vs. t were 
then carried out to find the damping time constant. This ap- 
proach was reasonably objective and robust, but can be criti- 
cized on the grounds that the core density is not constant and 
the density center is moving with time (Fig. |9]l, making the 
interpretation of E^no problematic. 

3. Given the difficulties with evaluating and interpreting the 
total energy of the BH, we chose in the end to quantify the en- 
ergy damping purely in terms of the the BH's kinetic energy. 
As noted above (Fig. [8] and associated text), V^{t) exhibits a 
nicely exponential decay with a well-defined time constant, 
and the decay is observed to continue over ~ 2 decades in ki- 
netic energy in the case of the simulation with the largest A^, 
until the BH's kinetic energy reaches the Brownian value. We 
evaluated the associated time constant by carrying out least- 
squares fits of InV^ to time, yielding the coefficients {Vf ,t) 
in the expression 

y^(f)wy/e"''"^'^/^. (20) 

Table |2] gives the r values derived from this method. We 
present results only from A^-body integrations with Vyck > 
0.4yesc since the smaller kicks did not excite distinct BH-core 
oscillations, as discussed above. To the extent that the mo- 
tion approximates a damped SHO, the energy damping time 
is identical to the time constant for decay of the kinetic en- 
ergy alone, and henceforth we will refer to r as the "energy 
damping time constant." However in practice, we will use 
equation (l20t only to predict changes in (V ^) . 

The energy damping times in Table |2] can immediately be 
scaled to physical units using equation (|2]i. Such a scaling pre- 
sumes that the core properties of our A'-body models - which 
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presumably determine t - are related to global properties in 
the same way as in real galaxies. A better scheme would re- 
late T directly to the parameters (pc, ctcMbh) that describe the 
conditions in the core. Since we do not understand the mech- 
anism(s) responsible for the orbital damping in Phase 11, we 
experimented with several ways of plotting r versus core pa- 
rameters. 

Figure [TO] shows that a reasonably tight correlation exists 

when ijJcT is plotted against al /(G^^^py^MBn)- This is the 
expected dependence if dynamical friction is responsible for 
the damping (cf. equation [TSll. However, the effective value 
of In A needed to produce the measured damping times is very 
small, 0.1 <lnA<0.3 (FigurefTOli. This is yet another way of 
stating that orbital decay in Phase 11 is much slower than pre- 
dicted by Chandrasekhar's formula - roughly a factor 10-20, 
if we adopt In A w 2.5 for the expected value of the Coulomb 
parameter (Fig. |5]l. 

In terms of this scaling. Figure [TO] allows us to express the 
damping times in Phase II as 

-'''>«V(5f^)""(5|j)' (22) 

where the second line uses the Mbh - cr relation 
jFerrarese & Fordll2005h . Based on Figure[8]and on the other 
arguments given above, we expect the scaling in equation ( |22] | 
to be independent of A^, i.e. of stellar mass. 
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Table|2]also gives estimates of 7}/, the elapsed time in Phase 
II. We defined 7}/ as the time, measured from the end of Phase 
I, required for the BH's velocity to fall to its rms value in 
Phase III, assuming the time dependence of equation ( l20b . Ta- 
ble|2]shows that Tu is typically longer than Tj. 

In a galaxy with ^10^ stars, V^^^^^ would be much lower, 
and Tii correspondingly longer, than in our models. Assuming 
that the exponential dependence of energy on time persists to 
arbitrarily low values of E, the additional time spent in Phase 
II would be 

T\n{N.^x/N) (23) 

where A^gai is the number of stars in the galaxy. We used 
the set of A^-body simulations in Figure [8] to test this depen- 
dence. According to equation (l23T l. doubling the number of 
particles should extend the elapsed time in Phase II by an ad- 
ditive amount of Tln2 = 0.693t w 2.2, given that the mean t 
value for the four integrations is 3.2. Figure |8] confirms this 
prediction for 0.25 x lO*" < < 2 x 10^ 

Accordingly, Table |2] also gives values of 7}/ calculated 
from this formula, for A^g^i = (3 x 10^3 x 10'",3 x 10",3 x 
10'^). Conversion from the A^-body units of Table |2] to years 
is discussed in §|6l 

The exponential nature of the damping implies that the dis- 
tribution of displacements during Phase II is approximately 
uniform in In Ar. 

4. EFFECTS ON THE STELLAR DISTRIBUTION 

The displacement of the BH due to gravitational radiation 
recoil affects the stellar distribution and therefore the den- 
sity profile of the host galaxy. We expect the stellar struc- 
ture inside the core to be particularly affected by the motion 
of the BH, with important implications for the shape of the 
brightness profile in the inner region. In order to evaluate the 
changes induced by the escaping BH, we constructed spatial 
and projected density profiles for all models at the end of the 
simulations, when the BH is well into the Brownian regime. 

Figure [TT| shows the space (left plots) and projected (right 
plots) density profiles in models Al, A2, B for Vyck = 
(0.1,0.3,0.5,0.9)Vesc- We constructed these density profiles 
using the kernel-based algorithm of Merritt et al. ( 2006a). 
Particle positions were first shifted to coordinates that placed 
the BH at the origin. The algorithm uses an angle-averaged 
Gaussian kernel and modifies the kernel width based on a pi- 
lot (nearest-neighbor) estimate of the density in order to main- 
tain a roughly constant ratio of bias to variance in the final 
density profile. The projected density Tj(R) was computed via 
numerical projection of the space density. In order to reduce 
the noise still further, we combined multiple snapshots at late 
times and performed the fit on the combined data sets. 

Figure [TT] shows that a large core develops in the simula- 
tions due to the escape of the BH and its several passages 
through the central region. As the BH oscillates under the ef- 
fect of the kick, it transfers energy to the surrounding stars, 
thus pushing them to larger distances. The stellar density 
in the core drops and the slope of the inner distribution de- 
creases, leaving an inner profile that is flatter than the initial 
one. The amount of flattening in the profile or, equivalently, 
the mass deficit with respect to the initial profile (shown in 
the figure with the black solid lines), increases monotonically 
with the kick velocity. 

It is interesting to assess whether the final A^-body profiles 
are consistent with the core-Sersic law, which is comm only 
fit to galaxies with evacuated cores dGraham et alj|2003l) . The 



TABLE 3 

Fit parameters for models Al, A2 and B. 



Vkick 


rb 


n 


a 


7 


Re 


Si, 


Al 


0.1 


0.013 


4.04 


3.1 


0.21 


0.93 


6.3 


0.2 


0.016 


4.05 


4.2 


0.24 


0.93 


5.6 


0.3 


0.017 


4.05 


4.1 


0.19 


0.93 


5.4 


0.4 


0.018 


4.06 


3.7 


0.14 


0.93 


5.1 


0.5 


0.019 


4.06 


3.6 


0.11 


0.93 


4.9 


0.6 


0.020 


4.06 


3.5 


0.08 


0.93 


4.6 


0.7 


0.021 


4.07 


3.1 


0.04 


0.92 


4.3 


0.8 


0.022 


4.07 


2.8 


0.02 


0.92 


4.1 


0.9 


0.024 


4.05 


2.9 


0.05 


0.93 


3.9 


A2 


0.1 


0.014 


4.04 


7.4 


0.34 


0.93 


6.5 


0.2 


0.015 


4.04 


6.7 


0.31 


0.93 


6.2 


0.3 


0.016 


4.04 


6.5 


0.25 


0.92 


5.9 


0.4 


0.015 


4.05 


3.5 


0.14 


0.92 


5.9 


0.5 


0.017 


4.05 


4.0 


0.12 


0.92 


5.4 


0.6 


0.019 


4.05 


4.2 


0.16 


0.92 


4.9 


0.7 


0.018 


4.07 


2.7 


0.08 


0.92 


4.9 


0.8 


0.022 


4.06 


3.6 


0.07 


0.92 


4.4 


0.9 


0.026 


4.07 


4.1 


0.11 


0.92 


3.9 


B 


0.1 


0.020 


4.05 


1.9 


0.16 


0.92 


4.4 


0.2 


0.026 


4.04 


2.9 


0.20 


0.93 


3.8 


0.3 


0.030 


4.04 


3.7 


0.20 


0.93 


3.4 


0.4 


0.034 


4.05 


4.3 


0.16 


0.92 


3.1 


0.5 


0.034 


4.06 


3.1 


0.12 


0.92 


3.0 


0.6 


0.035 


4.06 


3.0 


0.09 


0.92 


2.8 


0.7 


0.039 


4.08 


3.0 


0.07 


0.91 


2.6 


0.8 


0.042 


4.08 


2.9 


0.05 


0.91 


2.4 


0.9 


0.044 


4.09 


2.6 


0.02 


0.91 


2.3 



core-Sersic law is: 



S(7?) = S 



1 + 



m 



-A[(«°+C)/«°] 



l/„a 



I]'=S,2-^/"/('"°'^/*0"". 



(24) 
(25) 



where is the density at the break radius ri, and the other 
parameters are as in equations (|2|i and To carry out the fits 
in a manner as similar as possible to the procedure followed 
by observers, we counted the projected particle positions in 
bins equally spaced in log/?. The parameters (/?e,'"b,Q!,M,S ) 
were then varied until the summed residuals in /i = -2.5 log E 
were minimized. 

The best-fit parameters for models Al, A2 and B are listed 
in Table |3] Three of the best fits for model Al are shown in 
Figure [12] (lines) together with the projected density profiles 
computed from the A^-body data (points). 

It appears that the host galaxies to recoiling BHs are well 
represented by core-Sersic profiles. In particular, the fits 
show, once again, that the core tends to expand as the BH 
oscillates in and out of it, and that the final core size scales as 
ri, ^ A^BH V'lf^jj, with 0.3s;/9S0.6. In addition, the transition 
from the inner power law to the outer Sersic profile is rather 
sharp, with best-fit values of a in the range 2 < a < 7. The 
initial « = 4 de Vaucouleurs outer slope is not substantially 
modified by the BH. 

A flatte ning of the inner profile is also observed in the simu- 
lations of lBovlan -Kolchin et al. (2004), who follow the evolu- 
tion of a spherical stellar bulge with a recoiling central black 
hole using an A^-body tree code. They find that the density 
profile of the system evolves as a consequence of the gravi- 
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Fig. 1 1 . — Space (left) and projected (right) density profiles for models Al, A2 and B and different values of the kick velocity: Vyck = 0.1 Vcsc (green/dotted), 
Vkick = 0.3 Vcsc (red/dashed), Vkic^ = 0.5 Vise (cyan/dot-dashed), Vyck = 0.9 Vcsc (blue/long dashed). The black solid lines represent the initial profile, which is the 
same for each value of Vyck- 



tational radiation recoil and flattens substantially. A core of 
size equal to the BH sphere of influence forms on a relatively 
short time-scale, and remains even after several dynamical 
times. A flattening of the profile is observed for recoil veloc- 
ities smaller and larger than the central escape speed, though 
an additional flattening is present if the black hole returns to 
the core after the ejection. 

A measurable signature of a recoiling BH is the mass 
deficit, the net mass re moved from the central regions 
(iMilosavljevic et aT]l2002h . Mass deficits produced by recoil 
will add to the depletion caused by the pre-existing BH bi- 
nary, which ejects stars from the core during close encoun- 
ters. The deficit produced by the binary is proportional to the 
mass of the binary, with only a weak de pendence on the mass 
ratio and the initial density distribution (lMerrittl2006l) . There- 



fore, a binary BH can only produce a deficit Mdef ~ Mbh- 
This could explain the peak i n the distribution of observed 
mass deficits atMdef/MBH ~ 1 (Graham 20041 iFerrarese et"aD 
|2006|) . The tail of the distribution, however, extends to val- 
ues of Mdef/MfiH ^ 5. While such large valu es might be ex- 
plained as successive mergers (lMerritj|2006l) . a recoiling BH 
represents an interesting alternative. 

We evaluated the mass deficits in the final A^-body models 
by computing the difference in stellar mass, enclosed within a 
sphere of radius rs, between the initial and final space density 
profiles. Given the fact that the deficits depend rather sensi- 
tively on the value of rj, we computed Mdef as a function of 
Tj for a number of models and kicks. In all cases, Mdef first 
increases rapidly with rj and then flattens out to an approxi- 
mately constant value. Based on such tests, we concluded that 
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Fig. 12. — Projected density profiles for model Al computed from the 
A'-body data (points), compared with best-fitting core-Sersic models (lines), 
for three different values of the kick velocity (Vjiick = 0.2,0.4,0.8 Vesc). The 
insert shows a zoom into the central region. 
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Fig. 13. — Mass deficits, as defined in the text, for the different runs: Al 
(black), A2 (blue), B (red). Dashed lines show power-law fits. 

the most appropriate values to use for the computation of the 
mass deficits were as follows: r, = 0.05 for model Al, 0.04 
for model A2 and 0.1 for model B. The results for all three 
models are shown in Figure [T3] Also shown are least-squares 
fits to F = aX'', where Y = Mdef/MsH and X = Vkick/Vesc- The 
best-fit parameters are: 

Model Al :fl = 4.83, ^7=1.59 
Model A2:a = 5. 08, fe= 1.75 
ModelB :a = 4.31, ^7= 1.90 (26) 



The largest kicks result in mass deficits as large as 4-5Mbh, 
which is consistent with the largest observed deficits (iMerrittI 
12006'). Our definition of mass deficits as the difference in inte- 
grated mass between initial and final profiles implies that our 
estimates do not take into account any depletion prior to the 
kick. One should therefore add the contribution from the bi- 
nary evolution phase (Mjet ~ IA^bh) to our measured deficits 
before comparing with the observed values. 

The sensitivity of Mdef to r,, which presumably is a feature 
of real luminosity profiles as well, suggests that a more objec- 
tive way be found to measure mass deficits. 

We compare the projected density profiles obtained from 
the A^-body simulations to the brightness profiles of a sample 
of early-type galaxies in the Virgo cluster observed with the 
Advanced Came ra for Su rveys (ACS) on the Hubble Space 
Telescope (Ferra rese et al . 2006). In this study, the authors 
find that, while simple Sersic models generally provide a 
good representation of the global galaxy profiles, the bright- 
est galaxies require a power-law component within a charac- 
teristic break radius and are therefore best modeled with core- 
Sersic profiles. 

We select two representative galaxies in the sample and 
compare their surface density profiles with each of the 27 fi- 
nal profiles obtained from the simulations (9 values of the kick 
velocities 0.1 . . .0.9 for each of the 3 models Al, A2, B). 

The brightest Virgo galaxy, VCC1226 (M49, NGC 4472), 
has the largest value of mass deficit (Mdef/MeH ~ 4) and has 
a Sersic index that is not too far from our A^-body models, 
« - 5.91 On the other hand, VCC 731 (NGC 4365) has a 
relatively small core and a typical mass deficit of ^ IMbh 
(^erritt 2006). For each galaxy, we scale the A^-body profiles 
to have the same ri, and I](ri) as the galaxy itself. 

Figure[T4]shows that the brightness profiles of both galaxies 
can be reasonably well fit by (at least) one of the A^-body mod- 
els. In particular, the profile of VCC 1226 is well fit by mod- 
els with Vkick > 0.4yesc « 550kms-' while VCC 731 is well 
fit by models with Vkick~0.iyesc ~ llOkms"'. This indicates 
that observed brightness profiles, and even the largest cores, 
can be well reproduced by the gravitational recoil kicks. 

5. EVOLUTION TIMES IN REAL GALAXIES 

Given a galaxy's effective radius 7?e and total mass Mgai, 
equations (|2]i relate our A^-body units to physical units. We 
adopt the sca ling relations der ived from the ACS Virgo clus- 
ter survey of ICdte et alJ (l2004l) betwee n Re. and absolute blue 
magnitude Mb for early-type galaxies. iFerrarese et al.l ( l2006h 
found, for Virgo E galaxies fainter than Mg « -20.5, a mean 
relation 



logioTJe = 0.144-0.05 (Mb + 20) 



(27) 



where is in kpc. (Brighter galaxies obey a different relation 
and are considered separately below.) We relate Mg to galaxy 
mass using Gerhard et al.'s (2001) expression for the mass to 
light ratio in the blue band: 



logio 



1.17 + 0.67 log 



10 



10"Lp 



(28) 

Equation (l28l l was derived from dynamical modeling of galax- 
ies with Mb > -22.5 and represents an average for the matter 
within the effective radius, including dark matter if present. 

* Most of the blight galaxies in the ACS sample have n ^ 7. 
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Fig. 14. — Surface brightness profiles for the Virgo galaxies VCC 1226 
(top) and VCC 731 (bottom) from the ACS sample compared to the A'-body 
profiles obtained from the best fitting of the three models. The different lines 
correspond to the 9 different kicks Vitick/Vcsc = 0. 1 . . . 0.9. 

Combining these relations gives 



7?e ~ 1.2 kpc 



sal 



10"' M, 



0.075 



(29) 



The dependence of /?e on Mgai is weak, a consequence of 
the low slope of the R^-Mb relation. However we note that 
the scatter in this relation is large (e.g. iFerrarese et al.ll2006l 
Fig. 136). 

Some fiducial values, and their implied A^-body scalings 
(from equation|2]i, are: 



M„ai = 3 X lO^Me 



:3 X 10'"M 







M„ai = 3 X 10" Me 



/?e= 1.1 kpc 

[r] = 1.0xl0V [y] = 110kms- 
7?e = 1.3kpc 

[r]=4.1xl0V [V] = 315kms- 
7?e = l.Skpc 

[T] = 1.7 X lOV W] = 910kms- 



The trend of decreasing [T] with increasing Mgai reflects 
the well-knovvn hig her density of more massive galaxies 
dGraham et al.l 120031) . The central escape velocities in our 
models are 2.0 < Vesc ^ 2.2 in A^-body units, corresponding to 
Vesc ~ 2.1 X [V] « 2000kms-i when scaled to a 3 x 10"Mq 
galaxy. This agrees well with escape velocities of bright E- 
galaxies deri ved from more detailed modeling (e.g. Fig. 2, 
iMerritt et a l." 2004). 

All of the times listed in Table |2] can be scaled to phys- 
ical units using these relations. Figure [15] shows the result 



for the three fiducial values of Mg^i. We have used equa- 
tion ( |23] ) to correct the measured 7}/ values to different values 
of A^gai = Mgai/m^ assuming ni-^ = Mq; we also show, as con- 
servative lower limits, the 7}/ times obtained directly from the 
simulations. 

Figure [15] seem to suggest that return times depend discon- 
tinuously on Vkick, since Phase II does not appear to exist in 
our simulations when Vkick ^ 0.3Ve,sc. As discussed above, this 
might not be true in simulations with larger A^, or in real galax- 
ies. In any case, for Vidck 0.4yesc, return times are dominated 
by the time spent in Phase II ("BH-core oscillations"). 

The brightest galaxies, Mb < -21, appear to obey a dif- 
ferent scaling relation between Rf. and Mr than the rela- 
tion ( i29] l given above (IFerrarese et al.ll2()06h . Furthermore, 
these bright galaxies are typically fit by Sersic indices in the 
range 5 < « < 10, larger than the value n = 4 adopted here 
for the A^-body models. On the other hand, the brightest E 
galaxies often have resolved cores with well-determined sizes 
and densities (cf. §111). Furthermore, equation (l22l l gives the 
damping time t in Phase II in terms of core properties alone. 

We define the Phase II return times for these galaxies as the 
time for the BH's energy to decrease from 



1 2 2 



2 ^ 2 



the BH's energy when it first re-enters the core, to 



1 7 



3 



2 Mr 



(30) 



(31) 



the Brownian energy, assuming an energy damping time con- 
stant of r; for the latter we take equation (l22l) . This time is 



'15 



1.2 X 10 V 



(32) 



250km s- 



lO^M 



0PC 



"(. 



BH 



VlO'^Me 



Af = ln 



'An 



1 Mbi 

1 MBHMgal 



Mgai M, 



(33) 



with F K,2 the form factor defined above, and we have again 
assumed m,^ = M0. 

Figure [16] shows estimates of r and Mr for the six bright- 
est galaxies in the AC S Virgo sample e xcluding M87, which 
has an active nucleus (ICdte et al.l 12004*). Of course, this fig- 
ure is only meaningful under the assumption that the BHs in 
these galaxies have received large enough kicks to remove 
them completely from the core, i.e. Vkick ~ lO^kms"'. But 
if this did occur. Figure [16] suggests that return times would 
be of order 1 Gyr. Such a long time is comparable with the 
mean time between galaxy mergers in a dense environment 
like the Virgo cluster. Hence, a SMBH might never return 
fully to the center before another SMBH spirals in. 

6. OBSERVABLE CONSEQUENCES 

6. 1. Likelihood of Large Kicks 
Kicks large enough to remove SMBHs from cores, Vyck ^ 



0.4yesc, range from ^ 90kms for Mgai 



750kms-' for Mg^i = 3 x 10" M, 



to 



3 X 10^ Mq, to 
lOOOkms"' for 



Mgai = 3 X IO'^Mq, based on the fiducial scaUngs in §|5] The 
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Fig. 15. — Return times of kicked BHs. These are plots of the W-body 
values given in Table|2] scaled to physical units using equations ^2) and j29t . 
Lower (filled) symbols: 7}; Middle (open) symbols: Ti + Tn, with Tu taken 
directly from the simulations; Upper (filled) symbols: Ti + T//, with T// scaled 
to/Vgai using equation f23\. (a) Mg^i = 3 X 10' Mq; (b) Mg^i = 3 x IO'^Mq; 
(c)Mg,i = 3x 1O"M0. 




Fig. 16. — Estimates of the return time in Phase II for si iperm assive black 
holes in the six brightest Virgo galaxies, excluding M87 (Cote et al. 20(j^. 
This plot assumes that the SMBHs have received a kick large enough to re- 
move them from the core initially. Lower (ope n) symbols show the energy 
decay time constant in the core, t (equation |32) , while upper (filled) symbols 
show A/'t, where Af is the estimated number of time cons tant s required for 
the BH's velocity to decay to the Brownian value (equation |33). 



most propitious configuration for the kicks appears to be an 
equal-mass binary in which the individual spin vectors are 
oppositely aligned and oriented parallel to the orbital plane 
dCampanelh et: alj|2007at iGonzalez et al.ll2007al) . Assuming 
this most favorable orientation, and setting the spins to their 
maximal values, the maximum kick (oriented parallel to the 
binary angular momentum vector) is believed to scale with 
binary mass ratio q = M2/M1 < 1 as 



i 6 X 10 km s 



(34) 



(ICampanelli et alJl2007ch . Mass ratios as small as ^ « 0.2 



can therefore result in kicks > 1000 km s 



While the 



assumption of near-maximal spins is probably not an ex- 
treme one (e.g. Shapiro 2005; G ammie et alj ^004), ori- 
enting the BHs with their spins perpendicular to the or- 
bital angular momentum may seem od d, particularly in gas- 
rich galaxies dBogdanovic et al.l l2007l) . However there is 
considerable circumstantial evidence that SMBH spin axes 
bear no relation to the orientation s of the gas disks that 
surround them (iKinnev et alj 120001; iGalUmore et alJ 120061; 
iBorguet et alJ 120071) and this is presumably even more true 
with respect to the directions of infalling BHs in gas-free 
galaxies. If SMBH spins do orient parallel with orbital angu- 
lar momenta, the maximum kick is more modest and contains 
contributions from both the "mass asymmetry" (Mi y^Mj) and 
from the spins. The two kick components, both of which are 
parallel to the orbital plane, are believed to be approximately 
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independent and to scale roughly as 



spin ' 



'■V2 



qHl-q) 
(l + q)' ' 

q' 
(1+q)' 



■(a2-qai)\ 



(35) 
(36) 



Vi « V2 « 10"* km s~ ' and a denotes a dimensionless 
spin, -1 < g < 1 (CampaneUieLal. 2007c; Baker et al. 2007; 
iLousto & Zlochower 2007). Vmass peaks at ^ 200 km s~' for 
^ « 0.4 while Vspin peaks at ^ 600 km s"' for q Ri l,ai = 
-a2 = 1 . In this less-favorable configuration, kicks could re- 
move SMBHs only from the cores of low-to-moderate lumi- 
nosity ga laxies. Estimates of the kick velocity distribution 
(e.g. Schnittman & Buonannol l2007) are extremely uncertain 
since they depend on the unknown distributions of SMBH 
mass ratios, spins and spin orientations. In what follows, we 
will focus on the consequences of kicks that are large enough 
to remove SMBHs from galaxy cores and to excite the long- 
lived oscillations that we described above. 

6.2. Offset and Double Nuclei 

iLauer et aTl ( l2005b identified five galaxies in which the 
point of maximum surface brightness is displaced from the 
center of the isophotes defined by the galaxy on large scales. 
All are luminous, "core" galaxies. Contour plot s for two of 
the ga laxies, NGC 507 and 1374 (Figs. 17, 18 of iLauer et all 
l2005h . look strikingly similar to the "Phase 11" isodensity 
plots in Figure |9] Displacements are cited for NGC 507 
(0".06 w 19pc), NGC 1374 (0".02 w 2.1pc), and NGC 7619 
(0".04 « llpc), all of which are of order the core radii in 
these galaxies. The five galaxies with offset nuclei comprise 
12% of the Lauer et al. "core" galaxy sample; no offset nu- 
clei were found among the "power-law" (non-cored) galax- 
ies. Several of the offsets are close to the resolution limit, 
and some offsets might go unobserved due to projection, so it 
is likely that offset nuclei are quite common in "core" galax- 
ies. If the offsets are produced by oscillations like those in 
Figure |9] the SMBHs in these galaxies would be located on 
the opposite side of the galaxy photocenter from the point 
of peak brightness. Phase 11 oscillations can also produce a 
"double nucleus" morphology (e.g. Figure |9] frame 8) with 
the BH located at either the higher or secondary peak. This 
is a reas ona ble m odel for the double nucleus in NGC 4486B 
(ILauer et al.lll996i) . since the two peaks are closely matched 
in brightness and are offset by similar amounts (~ 6pc) from 
the galaxy photocenter. Galaxies with central mi nima in the 
surface brightness (e.g. NGC 4406, NGC 6876; ILauer et al.1 
(|2002)) might also be explained in this way. This model is 
probably not as appropriate for the more famous double nu- 
cleus in M31, since M31 is not a "core" galaxy, and one of 
the brightness peaks (the one as sociated with the S MBH) lies 
close to the galaxy photocenter (iLauer et al.lll993l) . 

6.3. Displaced AGN 

An ejected SMBH can appear as a spatially or kinematically 
displaced AGN (Kapoor 1976, 1983a b). A recoiling SMBH 
retains gas that is orbiting around it within a distance 



'"eff = 



GMbh 

•^kick 



!0.5pcM8yk; 



1000 



(37) 



with Mg = Mbh/IO'^Mq and Vk.iooo = Vkick/lOOOkms"!. An 
accretion disk if present would mostly be retained, and for 



kicks < lO-'kms"', rgff is large enough to encompass most 
of the broad emission-line region gas as well. Narrow emis- 
sion lines originate in gas moving in the gravitational potential 
of the host galaxy a nd woul d not fol l ow a r ecoiling SMBH 
(iMerritt et al.l l2006b). Bo nning et al.l (l2007h used this ai-gu- 
ment to search for kinematic offsets between spectral features 
associated with the broad- and narrow emission line regions. 
No convincing cases were found. This may be a consequence 
of the rapid decrease in SMBH energy during Phase 1 (Fig.O. 
In the longer-lived oscillations that characterize Phase II, the 
rms velocity of the SMBH drops from 



'90kms" 



103 M0 pc- 



1/2 



30pc 



when it first re-enters the core, to 



0.03 km s" 



Mr 



108 M,. 



-1/2 



200 km s- 



(38) 



(39) 



in the Brownian regime. Such small velocity offsets would 
be difficult to detect. An alternative approach would be to 
search for linear displacements A/? between the AGN emis- 
sion and the peak of the stellar surface brightness. This dis- 
placement is ^ fc at the start of Phase II, dropping to ^ 
y^m^^MBHrc in the Brownian regime; the exponential nature 
of the damping implies an approximately uniform distribu- 
tion of In A/? during Phase II. Relatively large (~ 10- lOOpc) 
offsets between the AGN and either the stellar density peak 
or the center of rotation have in fact been claimed in a 
number of galaxies ba sed o n integral-field spectro scopy (e.g 
iMediavilla & Arribaslll993l;lMediavilla et al 



i ectro sc 
200i. 



6.4. Wiggling Jets 

During Phase II, the SMBH oscillates sinusoidally within 
the core with roughly constant period. 



27r 

— w 1.4x lO'yr. , 

UJr V IO^Mq pC-3 



-1/2 



(40) 



and with velocities as given above. Such motion will induce 
periodic deviation s in the v elocity and direction of a jet emit- 
ted by the SMBH (|KaasEa& Roos 1992). Ifthejet is oriented 
perpendicularly to the direction of motion of the SMBH, the 
jet direction is fixed, and the jet material moves on a cylin- 
drical surface with radius equal to the radius of the SMBH's 
orbit. If the jet velocity has some component parallel to the 
SMBH's motion, the two velocities add a nd the cylinder be - 
comes a cone over which the jet precesses dRoos et al.l ll993'). 
Such models have been used to explain the heUcal distortions 
observed in a number radio sources; the inferred orbital peri- 
ods are typically 1 - lOOyr, and the jet accelerations are usu- 
ally ascribed to the orbit of the jet-producing SMBH around 
a second SMBH in a close (<C Ipc) binary pair. However 
some sources are fit by models with longer periods. For in- 
stance, the morphology of the C-type source 3C 449 has been 
reproduced assumin g jet forcing with a period of ~ 10^ yr 
dHardee et al.lll994l) . Such long periods are sometimes ex- 
plained in terms of bulk motion of the galaxy hosting the ra- 
dio source (Blandford & Icke 1978), but oscillations of the 
SMBH within the core might provide a tenable alternative in 
some cases. 

6.5. Oversized Cores and Hypermassive Black Holes 



Ejection of Supermassive Black Holes from Galaxy Cores 
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Cores generated by kicked SMBHs can be substantially 
larger t han those produced by "core scouring" from a binary 
SMBH jMilosavlievic & Merrittll200ltlMerri tl"2006'). partic- 
ularly when Vkick ^ O^Vesc- As shown in §|4] (Figures [TT} 
O, kick-induced cores can be as large as those observed in 
some of the brightest "core" galaxies, having mass deficits of 
4-5Mbh and core radii several times the SMBH influence 
radius, or ~ 5% of the galaxy's half- light radius. (Similar 
conclu sion s were reached alre ady by iBovlan-Kolchin et alj 
(|2P04) and iMerritt et all (12004 V) While the majority of ob- 
served mass deficits lie in the range 0.5 < Mdef/MBH S 1-5, 
some E galaxies have Mdef/ Msn ^, 3, to o large to be easily 
explained by core scouring. Lauer et aP (l2007h invoked the 
oversized cores, along with other circumstantial evidence, to 
argue that the SMBHs in the brightest E galaxies are "hyper- 
massive," Mbh }t 10'*'Mq. An alternative possibility is that 
the largest cores have been enlarge d by kicks. F igure[l3l com- 
bined with earlier A^-body results ( Merrittl2"00 6). suggests that 
the total mass deficit generated by a binary SMBH following 
a single galaxy merger is 

A^def = A^def.bin+A^def.kick 
Mdef,bin~0.7^"-2MBH, 
Mdef.kick ~ 5Mbh (Vkick/Vesc) ''^^ (41) 



where Mdef.bin and Mjef.kick are the mass deficits generated 
by "core scouring" and by the kick respectively and q = 
Mt/Mx < 1 is the binary mass ratio. It has been argued 
(Merritt 2006) that the ratio Mdef.bin /A^bh increases in mul- 
tiple mergers, and the same is likely to be true for kick- 
induced core growth. Thus, the decrease in typical values of 
of Vkick/Vesc with increasing galaxy luminosity might be off- 
set by the greater number of mergers that contribute to the 
growth of luminous galaxies, leading to comparable values 
of Mdef/MBH. In any case, the possibility that core growth is 
dominated by the kicks should be considered in future studies. 
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